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We present single and binary black hole simulations that follow the "moving puncture" paradigm 
of simulating black-hole spacetimes without excision, and use "moving boxes" mesh refinement. 
Focussing on binary black hole configurations where the simulations cover roughly two orbits, we 
address five major issues determining the quality of our results: numerical discretization error, 
finite extraction radius of the radiation signal, physical appropriateness of initial data, gauge choice 
and computational performance. We also compare results we have obtained with the BAM code 
described here with the independent LEAN code. 
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I. INTRODUCTION 

More than thirty years after the first numerical sim- 
ulations of binary black hole dynamics P, 0, the nu- 
merical relativity community is now ready to compare 
binary black hole simulations with experimental data. A 
series of recent breakthroughs [E 0> IE IE nas l e& d 
to a phase transition in the field: long-term evolutions 
of inspiraling black holes that last for one orbit and 
more have been obtained with several independent codes 
@SEHIIEIUIHm,and accurate gravitational- wave 
signals have been computed. 

Coincidentally, these breakthroughs parallel the first 
science runs of large-scale interferometric gravitational- 
wave observatories at design sensitivity The inspi- 
ral and merger of black-hole binaries are considered to be 
among the most promising sources for this current gener- 
ation of Earth-based gravitational wave detectors, and it 
has become feasible for numerical relativity to contribute 
to the analysis of experimental data. 

Such contributions will require large-scale parame- 
ter studies, and correspondingly large computational re- 
sources. The crucial current technical problem in the 
field is the efficiency of the simulations, and the estab- 
lishment of a "data analysis pipeline" , connecting analyt- 
ical calculations of the early inspiral and late ring-down 
phases with numerical simulations, and especially with 
gravitational-wave searches in actual detector data. 

In this paper we present a new version of the BAM 
code [E El for binary black hole simulations that 
follows the "moving puncture" paradigm of simulating 
black-hole spacetimes without excision, and use "moving 
boxes" mesh refinement. We give a detailed presentation 
of our methods, which will serve as a reference for future 
work, and we use simple test cases of single and orbiting 
black holes to calibrate our methods. 

We give a detailed discussion of convergence and ac- 
curacy of our code, and address further issues deter- 
mining the quality of our results: finite extraction ra- 
dius of the radiation signal, physical appropriateness 



of initial-data parameters, gauge choice and computa- 
tional performance. We compare evolutions from dif- 
ferent initial-data parameters and conclude that Post- 
Newtonian methods provide excellent estimates for initial 
positions, momenta and masses for quasi-circular orbits 
in the non-spinning equal-mass case, removing the ne- 
cessity of complex initial-data studies. We present new 
results concerning the damping parameter r\ in the pop- 
ular F-driver shift, which was originally introduced to 
handle long-term coordinate drifts, but is now found to 
also have the beneficial effect of magnifying the size of 
the apparent horizons and thus changing the spatial res- 
olution requirements. We also present timing and perfor- 
mance results: we have been able to perform some of our 
highest resolution runs at a computational cost of ss 500 
CPU-hours (see Table [HJ , giving rise to the hope that 
numerical relativity will indeed be capable of large-scale 
parameter studies. 

In Sec. [n] we recall the basic equations of the moving- 
puncture method, followed by a detailed description of 
our wave extraction algorithm and computation of ADM 
and Bondi quantities in Sec. 11111 Our numerical meth- 
ods and code structure are presented in Sec. IIVI In 
Sees. IVII and IVIII wc describe our results for single 
black holes, orbiting black holes evolved using standard 
quasi-circular-orbit initial-data parameters (allowing di- 
rect comparison with the LEAN code [HI), and orbiting 
black holes with alternative initial-data parameters. We 
conclude with a discussion in Sec. IVIIII 



II. THE PUNCTURE METHOD AND MOVING 
PUNCTURES 

A. Initial Data 

We will model A^-black-hole initial data by adopting 
the Brill-Lindquist wormhole topology ^3] with N + 1 
asymptotically flat ends for our initial geometry, thus 
enforcing the presence of N "throats" . The asymptoti- 
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cally flat ends are compactified and identified with points 
Ti on i? 3 . The coordinate singularities at the points 
resulting from compactification are referred to as punc- 
tures. In the context of initial data, punctures have been 
extensively studied and the treatment of the singularity 
in the constraint equations is well understood |18j . From 
a physical point of view the puncture representation of 
black-hole initial data is particularly appealing because 
it provides a simple prescription for associating masses, 
momenta and spins with any number of black holes. 

Initial data consist of the positive-definite metric and 
extrinsic curvature (gij, Kij) induced on a spatial hyper- 
surface S with timelike unit normal n l . We choose our 
sign conventions as n l rii = — 1, 

K a b = —■^—£n9ab- 

We construct these data using the conformal transverse- 
traceless decomposition of the initial- value equations, 
outlined in |l9l , and related to other conformal decompo- 
sitions in p0| . The spatial metric and intrinsic curvature 
are conformally related to counterparts on a background 
space via an initial conformal factor ipo, and the confor- 
mal extrinsic curvature is split into trace and trace-free 
parts: 



9ij 



-,9ij 



A, 



(1) 
(2) 



where K — g'^ and Ay is trace-free. 

We choose an initially flat background metric, gfy — 
Sij, and a maximal slice, K = 0. The second choice 
decouples the Hamiltonian and momentum constraints. 
The momentum constraint now takes the form 



djAij 



0. 



(3) 



and admits the Bowen-York solutions 12 lj for any number 
of black holes with prescribed ADM linear and angular 
momenta. 

The Hamiltonian constraint becomes an elliptic equa- 
tion for the conformal factor with a solution of the form 

mm mm mm 



IpO = '4>BL + U, 
4>BL = 1 



N 

E rrij 



(4) 
(5) 



The function u is determined by an elliptic equation on 
B? and is C°° everywhere except at the punctures, where 
it is C 2 . The parameterize the mass of each black hole, 
but actually equal the total mass of the black hole only in 
the special case of the Schwarzschild spacetime. In the 
general case, the ADM mass at the ith asymptotically 
flat end (i.e., the puncture) is given by 



Mi 




(6) 



where Ui is the value of u at the ith puncture, and dij 
is the coordinate distance between punctures i and j. 
This quantity has been found to agree within numerical 
uncertainty with the apparent-horizon mass Maha for 
non-spinning punctures |2jJ, and for spinning punctures 
we have found it to agree with the black-hole mass given 
by a modification of the Christodoulou formula [2F 
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AHA) 



(7) 



Throughout this paper, lower-case m s will refer to the 
mass parameter in the ansatz (JSJ, and M\ will refer to 
the black- hole mass given by Eq. ©. When we desire 
particular values of Mj, we make initial guesses of rrii by 
inverting Eq. JfJJl, and iteratively improve the rrn based 
on successive values of m until the Mi equal the desired 
values to within 0.02%. We will denote by M the to- 
tal black-hole mass of the system under investigation, 
and typically use M as the mass scale for quoting results 
(e.g. when to express time or distance in terms of a mass 
scale). 

To complete the definition of the initial data, we also 
need to specify initial values for our gauge quantities, the 
lapse function a and shift vector /3\ At time t — we 
use 



a = 1 or a = ip 
f3 l = 0. 



(8) 
(9) 



Both choices for the lapse have been used successfully, 
although the "pre-collapsed" lapse suggested in 0, |2jj is 
found to reduce initial gauge dynamics and is our stan- 
dard choice. Both lapse and shift are updated by evo- 
lution equations depending on the physical variables, as 
described below. 



B. Evolution System 

We evolve the initial data with the BSSN system [3(], 
l3ll ]. On the initial slice the standard BSSN variables <fr, 
gij, Aij, K, and P are related to the variables in the 
conformal transverse-traceless decomposition by 



(j) = InV'o 
A-ij = tp A^ 
r = -dfg l \ 



(10) 

(11) 

(12) 



and gij and K are unchanged. These variables are 
evolved using 



d (f> = ~a aK > 
6 

-a(KAij - 2A lk A k 3 ), 



doAij = e-^l-DiDja + aRij] 



TF 



(13) 
(14) 

(15) 
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d K = -D'Dia + a(A i:j A tj + \k 2 ), (16) 

-f 3 dj(3 % + -FdjPi - 2A ij d j a 

+2a (r) k ti k + eAVdrf - l^Kj ,(17) 

where da — dt — Cp, D% is the covariant derivative with 
respect to the conformal metric (jij, and "TF" denotes 
the trace-free part of the expression with respect to the 
physical metric, X7. F = X^ — ^gijX k . The Ricci tensor 
Rij is given by 

Rij — + Rfj (18) 

Rij = —^9 lm dldm9ii + 9k(idj)^ k + f fe f (ij)fe + 

~lm ( 2 ff (< f j)jm T f*mf Hi ) , (19) 

Rt = -2D i D j ct)-2g ij D k D k (t) + ADi(j>D j (t>- 

4g ij D k cf>D k ct>. (20) 

The Lie derivatives of the tensor densities 4>, and Ay 
(with weights 1/6, —2/3 and —2/3) are 

6 

2 

^/35ij = 9ijd k gij + gikdjfi k + gjkdif3 k - -hflkP , 

2 

CpAy = A,,ih A,j + A lk dj(3 k + A jk diP k - -A^d^. 

It is common to evolve the BSSN system as a par- 
tially constrained scheme, where one or both of the al- 
gebraic constraints det(g) — 1 and Tr(Aij) = are 
enforced at every full or intermediate time step of the 
evolution scheme. This has been found to be necessary 
to obtain stable, accurate evolutions of black-hole punc- 
tures in several cases, see e.g., [2^, [S^] . Likewise, the 
algebraic constraints and the first-order differential con- 
straint F' = —djg*i can be used for the source terms 
without affecting well-posedness, but changing for exam- 
ple the source terms of the constraint-propagation equa- 
tions. In the BAM code we make the following choices: 

• Wherever F l appears undifferentiated, we explicitly 
use r l = —djg l J instead of the evolved variable T\ 
Otherwise, F 1 is used. 

• The algebraic constraints Tr(Aij) = and 
det(gij) = 1 are imposed whenever the right-hand 
sides are calculated, and also at the end of each evo- 
lution step. (Imposing the algebraic constraints at 
each evolution mini-step does not imply that they 
will hold after each full time step, because of the 
nonlinearity of the expressions involved.) 



Further important choices concern the treatment of the 
conformal factor and the gauge choice. We will describe 
below the most popular choices that have been used suc- 
cessfully in the literature, and present some comparisons 
of different choices in Sections Mand IVII 

C. Choices for the conformal factor 

Let us first recall the fixed-puncture method, where the 
puncture pole is treated analytically and the punctures 
do not move. This is described in detail in 29]. The 
BSSN conformal factor is split according to 

(p = hxip = ^ + laipBL, (21) 

where £ is assumed to be regular at the puncture. In an 
evolution the regular function £ is evolved via the cor- 
responding version of Eq. (|13fl . and the logarithmically 
singular part In tpBL is kept constant. The key issue in 
the whole approach is to show that all evolved variables 
remain sufficiently regular at the punctures during evo- 
lution. In [2|j, a detailed analysis is given in terms of 
power counting arguments. 

The disadvantage of this method is that it assumes a 
natural split of the conformal factor according to Eq. 121|) 
throughout an evolution, i.e., that the slices remain con- 
nected to all asymptotically flat ends. 

In the moving-puncture approach the entire confor- 
mal factor is evolved. No assumption is made about the 
geometry of the slices. The slices are now allowed to 
approach whatever geometry is preferred by the gauge 
conditions. It turns out that in this preferred geome- 
try the conformal factor does not maintain the 1 jr Brill- 
Lindquist p ole, and instead develops a \ j\fr pole at the 
"puncture" [33J . The puncture ceases to represent a sec- 
ond infinity, and instead corresponds to a surface inside 
the horizon. The space outside this surface can be ac- 
curately resolved with a finite-difference code. We can 
then regard the moving-puncture method not as a mere 
trick to prolong the lifetime of a black-hole evolution, 
but rather as an elegant and simple alternative to ex- 
cision techniques: the singularity is not cut out of the 
numerical grid, it is avoided by the choice of gauge [33| . 

The question now is how to evolve the divergent confor- 
mal factor. In practice two proposals have been found to 
work, which we will call the ^-method and the ^-method. 
In the 0- method , one works directly with the original 
BSSN variable <f>, 

<f> = hxip, (22) 

and the evolution system remains as Eqs Q13 )) -Q17 [l . 
The purely experimental result is that finite differenc- 
ing across the ln(r) singularity at r = leads to stable 
evolutions. 

In the x _m ethod |3j , a new conformal factor is defined 
that is finite at the puncture, 

X = (23) 
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doX = ~x(aK - djF), (24) 

Now Eq. (|24|l replaces Eq. (|13fl in the evolution system. 
If ip has the usual 1/r pole at the puncture, then \ = 
0(r 4 ) at the puncture. As discussed in j^jj, the behavior 
changes to \ = 0(r 2 ) during the evolution. 

This approach does not rely on finite differencing of a 
singularity, but the singular structure of the black hole is 
incorporated in the vanishing of x a t the puncture. Be- 
cause of divisions by x present in the evolution equations, 
care needs to be taken in the numerical implementation 
to avoid divisions by zero or discontinuities arising out 
of unphysically negative values of x- We find that these 
problems can be avoided if x is consistently replaced in 
the right-hand sides of (|T31) - (fr7|) by max(x,e) (for some 
small e) wherever divisions by x occur. As a general rule, 
we choose e as follows. We know that near the puncture, 
X ~ (r/2m) 4 in the initial data, and later evolves to the 
form x ~ {r/ra) 2 . We therefore expect that x w iU n °t 
fall far below its initial minimum value, and choose e to 
be less than (r min /2m) 4 . 

D. Choices for the gauge 

The second ingredient in the moving-puncture method 
is a modification to the gauge choice. Both approaches 
now rely on the "covariant" form of "1+log" slicing [34| . 

(d t - (3 l d t )a = -2uK. (25) 

The shift advection term had been dropped in the ver- 
sion of "1+log" slicing used in the analytic fixed-puncture 
approach: 

d t a = -2aK, (26) 

and also in the first version of the moving-puncture ap- 
proach of . An attractive feature of Eq. l(2l)| is that the 
slicing is asymptotically maximal for a stationary solu- 
tion, such as the final Kerr black hole of a merged binary. 
However, Eq. (|26|l admits undesirable zero-speed modes 
in the BSSN system jsEHil and Eq. (23) turns out to be a 
better choice for moving-puncture evolutions The 
stationary Schwarzschild slicing with Eq. (|25() is given in 
[21E3. 

For the shift, we use a gamma-freezing condition. The 
original gamma-freezing condition introduced in |29| is 

d t P = d t B l = d t r - V B\ (27) 

Variants of this condition 0, ^J, consist of replacing 
some or all of the dt derivatives with do = dt — (3 l di. 
We will label these options with reference to each of the 
three time derivatives in J57J): "ttt" denotes that dt is 
used for all three derivatives, "Ott" denotes that <9 is 
used for the first time derivative, and d t for the other 
two, and so on. The properties of the different choices 



are studied in j3f| |3(|. Reference j3|J proves that the 
combination of the BSSN equations with the "1+log" 
slicing condition (|25|l and the 000 shift choice is strongly 
hyperbolic in the sense of first order in time, second order 
in space systems [H, EHI E3> EH an d thus yields 
a well-posed initial-value problem. For the final results 
presented in Section IVII we quote results obtained with 
the "ttt" and "000" options, which are both found to 
yield stable evolutions. All our recent work is based on 
the manifestly hyperbolic choice "000", i.e., we make the 
replacement dt — > do everywhere. 

III. ASYMPTOTICS 

A. Using ^4 for wave extraction 

Extracting physical information from numerical simu- 
lations in general relativity represents a highly non-trivial 
task for two reasons. First, most of the functions nu- 
merically computed in the course of an evolution are in- 
herently coordinate dependent. Second, quantities com- 
monly used for the description of local systems in other 
areas of physics, such as energy and angular momentum, 
are hard to define in an unambiguous way in correspond- 
ing scenarios in general relativity. For the problem at 
hand, the most important physical information to be ex- 
tracted are the energy and momenta radiated away in 
the form of gravitational waves and the precise shape of 
these gravitational waves as seen by a detector at large 
distances from the source. 

In the past, the extraction of these quantities from nu- 
merical simulations has been performed using either the 
Zerilli-Moncrief (see e. g. or the Newman-Penrose 
approach. In this work we focus on the calculation of 
the Newman-Penrose scalar ^4. This method has been 
discussed frequently in the literature, but we provide a 
detailed description to make clear the conventions we use 
(which can lead, for example, to differences in signs or 
constant factors of two), and to provide a complete, self- 
contained account of all the steps involved in calculating 
waveforms as well as radiated momenta and energy from 
the numerically evolved variables. For this purpose we 
will assume as known on a given hypersurface t = const 
the ADM variables </y and Ky. 

The Newman-Penrose scalar ^4 is defined by 

*4 = -R a ^n a m?n<m\ (28) 

where R a f3j8 represents the four-dimensional Riemann 
tensor (with the sign convention of an d n > ™ f orm 
part of a null-tetrad £, n, to, to. Specifically, £ and 
n denote ingoing and outgoing null-vectors whereas the 
complex- valued m is constructed out of two spatial vec- 
tors orthogonal to I and n, such that 

— £ ■ n — 1 = to • m, (29) 

and all other inner products between the four-vectors 
vanish. ^4 transforms as a spin-weight —2 field (that is, 
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under tetrad rotations which leave I and n unchanged 
and rotate m and ifi by an angle 9, we have ^4 — > 
e -2* e v[/ 4 ). Such objects represent symmetric trace-free 
tensor fields on a sphere (in our case R a p 1 $n a n 1 ) in terms 
of a complex scalar field. For a quick introduction to spin- 
weighted fields see e.g., |45j . There remains freedom in 
the choice of tetrad used in defining ^4. Here, we first 
construct a triad of orthonormal spatial vectors by ap- 
plying the Gram-Schmidt orthonormalization procedure 
to the three-dimensional vectors 



u l = [-y, %, 0] , 

v l = [x, y, z] , 
w* = g ia e abc u a v b . 

The tetrad vectors are then given by 



1 

V2a 



1 

y/2a 



1° 

m° = 



i 1 (-P 1 

V2\ a 

m l = —= (u l + iw l 
v2 



(30) 

(31) 
(32) 
(33) 



Next, we need to express Eq. <|28[) in terms of the three- 
dimensional quantities available on each time slice. This 
is achieved by virtue of the Gauss-Codazzi and the 
Mainardi equations which relate the space-time projec- 
tions of the four-dimensional Riemann tensor to its three- 
dimensional counterpart and the ADM variables accord- 
ing to 

_Li? a /3 7 <5 = IZafjyS + K al Kp6 — K a sKfj 7 , (34) 

-i-R f i/3"fSn fl = D^Kps - DgKpy, (35) 
IRtfriWn" = TLps-Kp^s+KKps, (36) 

where TZ a f3jS is the three-dimensional Riemann tensor, 
n a the timelike unit normal vector associated with the 
foliation and we follow York's notation for the pro- 
jection operator _L, which is, for example for an arbitrary 
tensor T a p, 



±-T af3 := (c^ Q + h»h a )(6 v + n v n p )T^. 



(37) 



In our coordinate basis adapted to the "3+1" decom- 
position, we are thus able to express \&4 exclusively in 
terms of the ADM variables as well as the triad vectors 
constructed from (|3"U|) according to 



* 4 



1 



(RabcdV a v c - 2±R abcd n a v c + LR abld n a n 1 ) 



(u b -iw b )(u d -iw d ). 



(38) 



The contributions of the individual modes i, m are ob- 
tained from projecting ^4 onto the spherical harmonics 
of spin weight —2. These projections are defined in 
terms of the scalar product 



A 



I III 



Jo Jo 



(39) 



which, in practice, is evaluated at a finite extraction ra- 
dius r ext . 

The spin-weighted spherical harmonics Y/ m can be de- 
fined in terms of the Wigner d- functions (e.g. |4(J) as 



Yf m (8,<p) = (-iy 

where 



2£ + l 

47T 



(40) 



%ns ifi) 

C 2 



(-1)* [{£ + m)\ (l - m)\ + s)\ 

t=Ci 



, n i/a 



(£ + m -t)\(t-s- t)\t\ (t + s- m)\ 

( n ic\\2l~\-m — s — It / • r\ /n\2t-\-s — m 

(cos (9/2) (sin 0/2) , 



(41) 



with C\ = max(0, m — s) and Ci = min(£ + m, t — a). For 
1=2 and spin- weight s — —2, we have 



Y" z - 

1 2-2 ~ 



y-2 _ 
r 2-l — 



1 20 — 



y-2 



y-2 
1 22 



64tt 



(l-cos^) 2 e- 2 ^, 



16tt 



svn.9 (1 — cos 9) e 



15 
327 



sm" 9, 



16tt 



sin^ (1 + cos6») e l 



64tt 



(l + cos<9) e 2 ^. 



(42) 



In practice, the integrand on the right-hand side of Eq. 
lUSH is evaluated on the Cartesian grid and interpolated 
onto a sphere of extraction radius r ext using fifth-order 
polynomials. The integration over the sphere is per- 
formed using the fourth-order Simpson method. 

While this procedure is straightforward from a numer- 
ical point of view, we emphasize one delicate point. In 
order to reduce the computational costs, numerical simu- 
lations are often performed with explicit use of symmetry 
properties of the spacetime under consideration. For this 
purpose it is important to take into account the trans- 
formation of the variables under inversion of the x, y 
or z coordinate. In our case the non-trivial operation is 
the symmetry across the xy plane. This problem man- 
ifests itself in the calculation of the modes according to 
Q39JI where the integrand is directly available only in the 
range < 9 < ir/2. Using the parity properties of the 
functions involved, however, we are able to transform the 
right hand side of Eq. P^jl into an integral restricted to 
the northern hemisphere. In particular, in the case of re- 
flection in the ^-direction ((x,y,z) — > (x, y, — z)), which 
is relevant here, the real part of "J 4 behaves like an even 
function, whereas the imaginary part of ^4 behaves as 
an odd function. 

Similarly, the harmonics Y^ 2 
complex conjugates of each other. In summary 



Y 2 _ 2 transform into the 



4-4(71- -0, 



4> 4 



(43) 
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Y ~ 



^22 V 



(44) 
(45) 



We use the following relation valid for arbitrary functions 
of 6 



f(6)d6 



tt/2 



\f(P) + f(*-O)]d0, 



IQ JO 

and are thus able to calculate 

r 2n r n/2 



* 4 F 22 2 sin 



(46) 



(47) 



(i 



^aY, 



4*22 



4*2-2 



Sill 



An equivalent change of basis to represent functions on 
the sphere has been discussed by Zlochower et al. in • 
In the study of numerical simulations of black-hole- 
binary systems, one is often interested in the amount 
of energy and momenta radiated away from the system 
in the form of gravitational radiation. In terms of the 
Newman- Penrose scalar ^4, these are given by the ex- 
pressions 



dE 
~dt 

dt 

dJz_ 
dt 





' r 2 r 




2 


lim 




1 ^ 4 dt 


dn 


r — >oo 


16^ J a 


J — 00 





lim 

r — ^00 




^ A dt 



0,1, 



■^idtdij dtt 



dCl 



ty 4 di 



(48) 
(49) 



where 



li = (— sin#cos</>, — sin 8 sine 



COS( 



(50) 



(51) 



We have listed these relations explicitly, because of dif- 
ferent conventions in use in the literature. In particular 
we emphasize the difference by a factor of 1/4 with Eqs. 
(22)- (24) of 0| which arises out of differences in the scal- 
ing of the tetrad vectors [cf. our Eqs. (|3*T|l - (f3*2")) with their 
Eq. (30)]. 

The expression for the energy can be simplified by us- 
ing the expansion of ^4 in modes £, m. Taking into ac- 
count the orthonormality of the spin-weighted harmonics 
we obtain 



dE 
~dt 



lim 



16tt 



E 

l,m 



Aprndt 



(52) 



In particular, this relation enables us to calculate the 
energy radiated in an individual mode. For the equal- 
mass systems considered in this work, we find > 99 % of 



the energy to be radiated in the form of the dominant 
I = 2, m = ±2 modes. 

Finally, let us note that in analyzing waveform modes 
as functions of time, it is extremely useful to split the 
complex function representing (or, say, the strain h) 
as 



rtf 4 (t) = A(i)e lv(t) , w(t) 



(53) 



as suggested in [lOj. In this paper this representation 
proves particularly useful to compare different initial 
data sets as in IVIII 



B. Total energy, linear and angular momentum 

In general relativity unambiguous notions of the 
energy-momentum four-vector and angular momentum 
can only be assigned to a spacetime as global quanti- 
ties, determined from the asymptotic structure of the 
spacetime. In this sense two types of quantity can be 
defined: those that are conserved by time evolution, and 
those that decrease with time, expressing the radiation 
of energy-momentum and angular momentum to infinity. 

The expression for the energy-momentum at spatial in- 
finity, which is time-independent and which corresponds 
to a four-vector under Lorenz transformations, was given 
first by Arnowitt, Deser and Misner in 1962 j49j in the 
context of the Hamiltonian formalism. This quantity 
is usually called the ADM energy-momentum, the time 
component being called the ADM energy or, somewhat 
inconsistently, the ADM mass, different from the rest 
mass to be defined below. The expressions can be given 
as limits of surface integrals defined at finite radius, and 
are evaluated in asymptotically Cartesian (regular) co- 
ordinates {x 1 } — where the components of the spatial 
metric tend to diag(l, 1,1) for large radii. The surfaces 
are then taken as spheres S r of radius r. 

We define the surface integrals (which we will also refer 
to as ADM integrals) 

E(r) = J s v^?V <Hj,k)dSi, (54) 



Pj{r) = ^J s V9(Ki-SiK)dSi 



(55) 



1 



Mr) = J ^ - KSl) dS t (56) 

which have to be evaluated in an asymptotically Carte- 
sian coordinate system. 

The ADM energy Madm and linear and angular mo- 
mentum Pj and Jj are then given by [l|| OHJ, yDj 



M ADM = lim E(r), 

r— »oo 

Pj = lim P,(r), 

r^oo 

Jj = lim Jj( r ) 



(57) 
(58) 
(59) 
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and the rest mass Mr can be defined as M R = M\ DM — 

For radiation processes we also require definitions of to- 
tal energy, linear and angular momentum that decrease 
as energy and linear as well as angular momentum are 
radiated to infinity. The appropriate quantities are the 
Bondi quantities [52 ] , which can be defined as taking the 
limit of the ADM inte gra ls not toward spatial, but rather 
null, infinity [5^, l55j|. i.e., the limit to infinite dis- 
tance is taken for constant retarded time instead of on a 
fixed Cauchy slice. In the context of our numerical treat- 
ment, the ADM and Bondi quantities can be computed 
rather accurately by computing values at several radii, 
and then performing a Richardson extrapolation (in ex- 
traction radius, not, as is more usual, in grid spacing) 
to infinity. Here the Bondi quantities can be computed 
at any time for a fixed extraction radius, and have to be 
compared between different radii by taking into account 
the light travel time between the timelike cylinders of 
different radii. This time delay can be estimated from 
a corresponding Schwarzschild solution as is done in j5|| 
by the difference in the values of the radial "tortoise co- 
ordinate" values as 

AT(R 1 ,R 2 ) = [R+ 2M\n{R/2M - 1)] B r Zr\ , (60) 

where the radii Ri are understood as Schwarzschild ra- 
dius (i.e. luminosity distance), and the Schwarzschild ra- 
dius can be estimated from the simulation's radial coor- 
dinate r by assuming it corresponds to the isotropic ra- 
dial coordinate in Schwarzschild spacetime, which yields 
R = r(l + 2M/r) 2 . 

IV. NUMERICAL METHOD 

The numerical method of our black-hole simulations 
is based on a method of lines approach using finite dif- 
ferencing in space and explicit Runge-Kutta (RK) time 
stepping. For efficiency, Berger-Oliger type adaptive 
mesh refinement (AMR) is used :57j. The new numer- 
ical results discussed in this paper were obtained with 
the BAM code [(j, Il5l llq , which implements a particular 
AMR strategy that we describe below (we also compare 
with published results obtained with Sperhake's LEAN 
code |K|). Although BAM also includes an experimen- 
tal oct-tree cell based algorithm that allows arbitrarily 
shaped refinement levels, this has not been used since a 
simpler box based algorithm is sufficient for black-holes 
binaries. 

The numerical domain is represented by a hierarchy of 
nested Cartesian grids. The hierarchy consists of L levels 
of refinement indexed by Z = 0, . . . , L — 1. A refinement 
level consists of one or more Cartesian grids with con- 
stant grid-spacing hi on level I. A refinement factor of 
two is used such that hi = h /2 l . The grids are properly 
nested in that the coordinate extent of any grid at level 
I, I > 0, is completely covered by the grids at level I — 1. 
Of special interest are the resolutions h ma x — Ziq of the 



coarsest, outermost level, and h m i n — h^-i of the finest 
level. 

Since we focus on the case of one or two black holes, 
a particularly simple grid structure is possible where 
each refinement level consists of exactly one or two non- 
overlapping grids. While the size of these grids could be 
determined by truncation error estimates or some field 
variable that indicates the need for refinement, for the 
purpose of convergence studies we have found it conve- 
nient to specify the size of the grids in advance. This 
allows, for example, the doubling of resolution within a 
predetermined coordinate range. Concretely, let Ni be 
the number of points in any one direction for a cubical 
box with Nj* points on level Z. On level I, center such a 
box on each of the black-hole punctures. If there are two 
punctures and the two boxes do not overlap, this is the 
layout that is used. If two boxes overlap, replace them by 
their bounding box, which is the smallest rectangular (in 
general non-cubical) box that contains the two original 
boxes. 

Assuming Ni = N (a constant independent of Z), a 
typical configuration around two punctures consists of 
two separate cubical boxes at I = L— 1, but for decreasing 
Z and increasing hi the size of the boxes increases until 
starting at some intermediate level the boxes overlap and 
a single rectangular box is formed, which towards Z = 
becomes more and more cubical. 

The hierarchy of boxes evolves as the punctures move. 
We use the shift to track the position x punc of a puncture 
by integrating 

^tXpmic ~ ~P ( x punc)> (61) 

cf. using the ICN method. The outermost box on 
level and also several of the next finer levels are chosen 
to be single cubes of fixed size centered on the origin to 
avoid unnecessary grid motion. 

Note that as long as one neglects the propagation of 
gravitational waves, the nesting described above repre- 
sents in a natural manner the 1 jr fall-off of the metric 
for a single puncture. For a single puncture and fixed 
N, doubling the grid-spacing going from level Z to I — 1, 
i.e., hi-i = 2hi, puts the boundary of a centered cube 
twice as far away. If a resolution of hi is sufficient to 
resolve the metric at 1/r, then 2hi should be sufficient 
to resolve the metric at 1/ (2r) since this is the slowest 
fall-off of any metric variable. This was the rationale for 
the nested box fixed mesh refinement (FMR) introduced 
in 0], which was found to work well in practice for the 
first 3d mesh-refinement evolutions of black holes [l6| . 
This FMR nesting strategy generalizes straightforwardly 
to the case of two moving punctures as outlined above. 

In the presence of gravitational waves further demands 
for spatial resolution arise: the wave amplitude falls off 
with 1/r, corresponding to the roughly constant ampli- 
tude of the "predicted" signal rip4 , while the wavelength 
is approximately constant. The spatial profile of the sig- 
nal thus requires constant radial resolution with increas- 
ing distance, while the amplitude fall-off leads to increas- 
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ing accuracy requirements as distance increases, in order 
to separate the waves from the background. Correspond- 
ingly, the grids need to be adapted when waves need to 
be traced accurately to typical wave extraction distances, 
which in a setup as presented here are still rather limited 
by computational cost. In actual runs it is thus conve- 
nient to use at least two different values for the Ni, one 
for the cubes that resolve the neighborhood of the punc- 
tures, another one for the levels where the wave extrac- 
tion is performed. For Berger-Oliger time-stepping most 
of the computational work is performed on the finest lev- 
els, so one chooses Nl_i as small as possible (while still 
covering the entire black hole with sufficiently fine reso- 
lution), and we can gain some extra resolution for wave 
extraction at small extra cost by using a larger box for 
the levels on which waves are extracted. 

The grids are cell-centered. For example, in one di- 
mension for the cell given by the interval [0, ho], the data 
on level is located at the point ho/2, on level 1 at 
h / 4: and 3/io/4, on level 2 at h /8, 3ho/8, bho/8, and 
7h /8, and so forth. Data is transferred between lev- 
els by sixth-order polynomial interpolation, where the 
three-dimensional interpolant is obtained by successive 
one-dimensional interpolations. 

On any given box with resolution hi, we implement 
fourth-order finite differencing for the spatial derivatives 
of the Einstein equations. Standard centered stencils are 
used for all first and second-order derivatives except for 
advection derivatives, /3 z <9j. For second-order finite differ- 
encing, the advection terms required one-sided differenc- 
ing for stability. For fourth-order finite differencing, we 
found that both centered and one-sided differencing can 
lead to severe instabilities with ICN time stepping, while 
"lop-sided" stencils lead to stable evolutions (cf. |58|). 
Our runs are performed using such lop-sided advection 
derivatives with fourth-order Runge-Kutta (RK4). 

The code allows us to add artificial dissipation terms 
to the right-hand-sides of the time evolution equations, 
schematically written as 

d t u -> d t u + Qu, (62) 

in particular we use the standard Kreiss-Oliger dissipa- 
tion operator (Q) of order 2r 

Q = <r(-h) 2r - 1 (D + y P (D_y/2 2r , (63) 

for a (2r — 2)-accurate scheme, with a a parameter reg- 
ulating the strength of the dissipation, and p a weight 
function that we currently set to unity. Adding artifical 
dissipation is apparently not required for stability in our 
runs, but we have used dissipation for RK4 evolutions to 
avoid high frequency noise from mesh-refinement bound- 
aries. We find that the inherently stronger dissipation 
of the ICN algorithm also rather efficiently suppresses 
noise from refinement boundaries, and our ICN test runs 
suggest that in this case the adding of Kreiss-Oliger dis- 
sipation is superfluous. 

All AMR results for two punctures reported so far are 
based on codes that involve at least some second-order 



component, while BAM in principle allows fully fourth- 
order AMR. In particular, we apply sixth order poly- 
nomial interpolation in space between different refine- 
ment levels so that all spatial operations of the AMR 
method are at least fourth order. However, there are 
three sources of second-order errors. One is the initial 
data solver, although this initial error appears to be neg- 
ligible in the cases we consider here. Another source of 
second-order error is the implementation of the radiative 
boundary condition. However, the nested boxes position 
the outer boundary at sufficiently large distances such 
that these errors do not contribute significantly (ideally 
because they are causally disconnected from the wave ex- 
traction zone). The final source of second-order error in 
our current runs is due to interpolation in time within 
the Berger-Oliger time-stepping scheme, which is worth 
discussing in some more detail. 

Berger-Oliger time-stepping can be stated as recursive 
pseudo-code for example as: 

evolveJiierarchy(Z, At) 
evolve(Z, At) 
if (I + 1< L) 

evolveJiierarchy (Z + 1, At/ 2) 

evolveJiierarchy (Z + 1, At/ 2) 
if (I > 0) 

restrict_prolong(Z — 1, I) 
regrid(Z) 

The recursion is started by calling the function "evolve_- 
hierarchy" for I — 0, i.e., beginning with the coarsest 
level. The function "evolveJiierarchy" evolves all levels 
from I to L — 1, the finest level, by a time step of At 
forward in time. First, level I is evolved by At by the 
function "evolve" . Then the function "evolveJiierarchy" 
calls itself recursively to advance level / + 1 and all its 
sublevels twice by At/2. The recursion ends if level I + 1 
does not exist, i.e., if I + 1 is not less than L, then 
"evolveJiierarchy" does not call itself again. Once all 
levels I through L — 1 have reached the next time level, 
information is exchanged between levels I — 1 and I, de- 
noted by a call to "restrict_prolong" if I > 0. In par- 
ticular, the refinement boundary of I is populated using 
information from I — 1. The result is the new level I. Fi- 
nally, the refinement hierarchy is updated by the function 
"regrid" . 

Although the time stepping used for evolution is 
fourth-order Runge-Kutta, there arises the additional is- 
sue of how to provide boundary values for the interme- 
diate time-levels of the Berger-Oliger algorithm that are 
not aligned in time with a coarser level (otherwise spatial 
interpolation can be used) . There are several options for 
fourth-order boundaries. 

The original suggestion by Berger and Oligcr is to in- 
terpolate in time (over several coarse levels at different 
instances of time) in order to obtain boundary values for 
a fine level. One can use three time levels of the coarser 
level to perform quadratic interpolation (third order in 
the time step) resulting in overall second-order conver- 
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gence when using a leapfrog scheme, e.g., as done in [l5j . 
However, the convergence order and the stability of the 
algorithm depends on the form chosen for the Einstein 
equations and on the time-stepping algorithm used. For 
example, quadratic interpolation for ICN and a first or- 
der in time, second order in space formulation can lead to 
a drop of convergence order and instabilities, see Schnet- 
ter et al. |(3l| . Other authors report success with different 
variants of time interpolation, e.g., |62l l63j| . 

An alternative approach is to replace the single point 
refinement boundary by a buffer zone consisting of sev- 
eral points, e.g., [(H], |64], Ell . The buffer zone approach 
can be expected to perform well for the transmission of 
waves through refinement boundaries, see e.g., (note 
that special methods like (6j| seem to achieve similar 
performance). The optimal number of buffer points is 
method dependent. For example, RK4 requires 4 source 
evaluations, and if the lop-sided stencil with 3 points in 
one direction is used, then the numerical domain of de- 
pendence for a given point has a radius of 12 points. 
Therefore, it is possible to provide 12 buffer points at 
the refinement boundary and to perform one RK4 time 
step with size 3 stencils that does not require any bound- 
ary updates. Only after the time step is completed, the 
buffer zones have to be repopulated. In the context of 
Berger-Oliger AMR, the buffer update is based on inter- 
polation from the coarser levels. Since every second time 
step at level I coincides in time with level Z — 1, one can 
provide 24 buffer points, perform two time steps, and 
then update the buffer by interpolation in space. With 
12 buffer points, one can interpolate in time to obtain 
data for the buffer points at intermediate time levels. 

For the simulations reported here, our standard setup 
is to use RK4 with dissipation and lop-sided advection 
stencils, 6 buffer points, quadratic interpolation in time, 
and Berger-Oliger time-stepping on all but the outermost 
grids. Let us comment on these choices. 

For some grid configurations we have encountered in- 
stabilities for very large, coarse grids, that experimen- 
tally are connected to the large time steps on the coarse 
grid. We were able to cure these instabilities by turning- 
off Berger-Oliger time-stepping for the outermost grids 
(cf. where this idea was introduced in a different con- 
text). 

To use fewer than 12 buffer points, we can interpolate 
into all buffer points before starting a RK4 update as 
described, and then evolve all points except the outer- 
most points located exactly on the boundary, which are 
kept fixed at their initial interpolated value. The inner 
points next to the boundary are updated using second 
order finite differencing for the centered derivatives and 
shifted advection stencils for the advection derivatives. 
Experimentally, using just 6 buffer points leads to very 
small differences compared to 12 buffer points, however 
smaller buffer zones lead to noticeable differences. Even 
though for large grids the number of buffer zones be- 
comes negligible, for the grid sizes that we have to use, 
the buffer points impact the size of the grids significantly. 



For example, for a box of size 64 in one direction, adding 
6 points on both sides instead of 12 points corresponds 
to a savings of 35% in the total number of points. For 
clarity, we always quote grid sizes without buffer points, 
because this is the number of points owned by a partic- 
ular grid. 

Using quadratic interpolation in time is, apart from 
the outer boundary treatment, the only source of second- 
order errors in the evolution scheme. We checked for 
a few cases that running without Berger-Oliger time- 
stepping entirely led to only small differences compared 
to other error sources. However, for sufficiently high res- 
olutions, quadratic interpolation in time should become 
the dominant error. In principle, we can resolve this is- 
sue by either not using Berger-Oliger time-stepping or 
by using larger buffer zones, which at the moment is pro- 
hibitively expensive in resources. 

We have also experimented with higher order in time 
interpolation, although a systematic analytical and nu- 
merical analysis beyond these first experiments is needed. 
Simply using additional coarse time levels was not suc- 
cessful. In general, if at time t a fine level I is not aligned 
in time with the coarser level 1—1, we use the grid func- 
tions on level I — 1 at different times to interpolate to 
time t. For quadratic interpolation these different times 
are t + At, t — At, and t — 3 At, where 2 At is the timestep 
on level I — 1. As mentioned before this kind of interpo- 
lation leads to overall second-order accuracy in time at 
the interpolated points. We routinely use this approach 
and it leads to stable evolutions. In order to obtain a 
fully fourth-order scheme we have included additional 
coarse levels at times t — 5 At and t — 7 At. However, this 
extended interpolation scheme over five different times 
leads to oscillations at the refinement boundaries, which 
are the points where we use interpolation in time. These 
oscillations increase with resolution and are thus likely 
instabilities which would cause the code to fail at suf- 
ficiently high resolution. At the resolution considered 
in this paper these instabilities do not cause the code 
to fail. However, they are a significant source of noise, 
which propagates out of the refinement boundaries into 
the rest of the grid. Since this noise is not convergent, 
it eventually spoils convergence in the entire grid. One 
reason for this problem may be the high degree of asym- 
metry in the interpolation stencil which uses four points 
before time t and only one after t. 

Finally, we note that BAM is MPI parallelized. The 
dynamic grid hierarchy with moving and varying boxes 
introduces an additional communication overhead com- 
pared to the FMR runs that BAM was used for previ- 
ously Q. For up to 128 processors scaling seems reason- 
able for a constant problem size per processor, but we 
do expect issues for larger processor numbers, which we 
have not been able to test yet. 
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V. SINGLE PUNCTURE WITH DYNAMIC 
CONFORMAL FACTOR 

A. Numerical experiments for a single stationary 
puncture 

In this section we apply the <p and x moving-puncture 
methods to evolutions of a single Schwarzschild punc- 
ture. This provides an excellent test case, because we 
can compare with the analytic results in 33], and study 
the convergence properties of the code without the added 
complication of moving mesh-refinement boxes. 

The initial data are as described in Section lil Al where 
U = and Aij = on the initial slice, and we choose 
mi = M = 1. We use a " pre-collapsed" lapse of a = i>^ 2 
for these runs, but stress that similar convergence proper- 
ties are found with an initial lapse of a = 1. The conver- 
gence series consist of evolutions with N — 64 3 , 96 3 , 128 3 
grid points in each box, and seven levels of refinement 
below the coarsest level (making a total of eight lev- 
els). The resolutions on the coarsest levels are h max — 
6M, 4M, 3M, and the resolutions on the finest levels and 
at the puncture are h min = 6M/128, 4M/128, 3M/128. 
The gauge choice is ttt, and r\ = 2.0/M. For these runs 
(and only for these) a uniform time-step was used on all 
levels (i.e., not Berger-Oliger) in order to fully test the 
fourth-order accu racy of the code. 

As discussed in [33j |. a 1+log evolution of Schwarzschild 
reaches a stationary slice, and at the puncture 1 = 
($if3 l = 0.5239 and the Schwarzschild radial coordinate 
is R = 1.3124M. After 50M of evolution, these values 
are reached to within 1.3% and 0.5% in the highest reso- 
lution runs using the method. With the \ method, the 
errors in /3 2 and R are 0.6% and 0.2%. Figure [I] shows 
several of the BSSN variables after 50M of evolution with 
the x method. 

The convergence of the 4> method is demonstrated in 
Figure which shows convergence plots of g xx , 4>, a, (3 V , 
and the Hamiltonian and momentum constraints. The 
data are taken along the y-axis at t = 50M on the finest 
level of the mesh-refinement scheme. The errors in the 
Hamiltonian and momentum constraints cover a wide 
range, so the logarithm of the scaled errors is shown. 
The differences are scaled assuming fourth-order conver- 
gence, and the code demonstrates good fourth-order con- 
vergence everywhere except at the points closest to the 
puncture. 

The <p variable shows extremely poor convergence at 
the puncture, but this is to be expected: <p diverges like 
ln(r) near the puncture. What is remarkable is that this 
non-convergent behavior remains localized at the punc- 
ture, and does not affect the accuracy or stability of the 
evolution as a whole. 

Figure [3] shows similar convergence plots for the x 
method. In this case the x variable, which should be- 
have like r 2 near the puncture, is seen to converge ev- 
erywhere. The constraints and (3 V also show better con- 
vergence properties near the puncture. This is consistent 
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FIG. 1: The BSSN variables g xx , x, a, and j3 y after 50M of 
evolution of a Schwarzschild puncture using the x method. A 
small pulse due to the initial adjustment of the gauge can be 
seen at about y = 60M in g xx . The main features of the other 
variables are confined to y < 10 M. 



with the comparison with the stationary 1+log solution, 
where we see that the x method was more accurate at 
the puncture. 

We draw three conclusions from these results. (1) Our 
code is fourth-order accurate for the resolutions used in 
this work, at least when the mcsh-rcfincmcnt boxes do 
not move and a uniform time-step is used. (2) The 
moving-puncture method extremely accurately reaches 
the stationary 1+log slicing, and, since the puncture no 
longer represents a second infinity, the solution is well- 
resolved up to the puncture. (3) Both the (f> and x meth- 
ods are stable and accurate, but the x method shows (as 
expected) better convergence properties at the puncture. 
As a test of the stability of the method at extremely high 
resolutions, we have also evolved a Schwarzschild punc- 
ture with resolutions of up to M/512 at the puncture, and 
found no signs of instability after 100M of evolution. 

In addition, we emphasize that the only variable that 
diverges at the puncture is </>. When the x method is 
used, all variables are finite at the puncture. Some vari- 
ables are discontinuous at the puncture. This leads to 
incorrect evaluation of finite-difference derivatives at the 
grid points closest to the puncture (the number of points 
depends on the width of the finite-difference stencil used), 
but these errors do not seem to propagate away from 
the puncture, and spoil the convergence of the variables 
in question only near the puncture. These errors could 
presumably be reduced or removed by using appropri- 
ate one-sided derivatives next to the puncture, but we 
have obtained sufficiently accurate results without need 
of such a sophisticated treatment. 
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FIG. 2: Fourth-order convergence of a Schwarzschild punc- 
ture after 50M of evolution, using the <j> method. Results 
were taken from runs with TV = 64 3 , 96 3 and 128 3 points with 
octant symmetry. The plots show the differences between the 
three runs, scaled to be consistent with fourth-order accuracy. 
For the Hamiltonian constraint and {/-component of the mo- 
mentum constraint, which should converge to zero, we show 
the logarithm of the scaled values. 
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FIG. 3: Fourth-order convergence of a Schwarzschild puncture 
after 50M of evolution, using the x method. The parameters 
match those used for the runs discussed in Figure [5] 
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FIG. 4: Left: coordinate location of R = 2M after 50M of 
evolution, as a function of the damping parameter, 77, with 
initial lapse a = 1 and a = ip~ 2 . Right: coordinate location 
R — 2M as a function of time, for Mr] = 0, 2.0, 3.5, using 
initial a = 1. 



term drifts in the metric variables. The effect of 77 in 
our new evolutions is demonstrated in Figure 01 which 
shows the coordinate location of the Schwarzschild hori- 
zon R — 2M after 50M of evolution, as a function of 77. 
We see that the coordinate size of the black hole differs by 
more than a factor of two between 77 = and 77 = 3/M; 
similar effects were alluded to in 0. As a result, differ- 
ent choices of 77 correspond to different effective numeri- 
cal resolutions across the black hole. For example, with 
7] = and a central resolution of M/16, there are about 
26 grid points across the interior of the black hole. With 
?/ = 3/M, there are about 59 grid points across the black 
hole — it is resolved twice as well! On the other hand, if 
the finest box in the mcsh-rcfinemcnt structure contains 
32 3 points, then this box contains the entire black hole 
when 77 = 0, but does not when 77 > 1.0/M. 

In any black-hole simulation, one must decide which is 
more important, the effective finest resolution, or the ef- 
fective size of the finest box. Perhaps more importantly, 
the effect of 77 on the coordinates shows that one must 
be careful when comparing runs that use different reso- 
lutions and/or box sizes, and different values of 77. 

Larger values of 77 also cause a larger drift in the hori- 
zon location with time. Although the geometry becomes 
stationary after about 40M of evolution, the numerical 
coordinates may not. This is clear from the lower plot in 
Figure^] We see that, if we wish to minimize the drift in 
the numerical coordinates, lower values of 77 are better. 
We will see similar results in Section IVII in the case of 
black-hole binaries. 



Numerical experiments for a single spinning 
puncture 



1. Coordinate dependence on n 

The geometry of the stationary 1+log slice is unique, 
but the coordinates of that final slice are not. One 
quantity that alters the final coordinates is the gamma- 
freezing damping parameter, 77. The parameter 77 was 
originally introduced in 29] for fixed-puncture evolutions 
to prevent oscillations in the shift vector as well as long 



We now look at results for evolutions of a single spin- 
ning puncture. These allow us to test the moving- 
puncture method for spinning black holes, and provide 
a non-trivial test of the wave-extraction algorithm for 
a black-hole spacetime. The initial data are now based 
on the Bowen-York extrinsic curvature for a single black 
hole with non-zero spin, which can be considered as a 
Kerr black hole plus Brill- wave radiation 21, 66, 6Hl68|. 
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FIG. 5: Real part of the I — 2, m = mode of r^4, extracted 
at r = 30mi. 



In an evolution the additional radiation will leave the 
system, and only the Kerr black hole will remain. The 
energy of the radiation has been estimated in the past 
by studying the initial data |66l l67| . with a radiation 
content of up to 3% for a near-maximally spinning Kerr 
black hole. 

We considered a Bowen-York puncture with mass pa- 
rameter mi = 1 and angular momentum parameter 
S z = 0.2mj. As discussed in Section fll Al the mass of 
the black hole can be estimated using Eq. For these 
data, the black-hole mass is M — 1.0155mj.. The Kerr 
parameter can then be estimated as a = s/M 2 — 0.194. 

The spinning puncture was evolved for IOOtoi using 
the <j) an d X methods. Convergence tests consisted of 
runs with seven levels of refinement, box sizes of 40 3 , 48 3 
and 64 3 points, and resolutions of the coarsest box of 
h m ax = 6mi, 5toi, 3.75mi. 

Figure[S]shows convergence plots for the <fi and x meth- 
ods. We have plotted the differences in Re(r& 4)20 be- 
tween the three grid sizes, and scaled the medium- fine 
difference by a factor of 1.57, consistent with fourth-order 
convergence. Both methods show reasonable fourth- 
order convergence for the first ~ 4QM of evolution, 
demonstrating that the wave-extraction algorithm is 
fourth-order convergent. Convergence in the waveform 
(and the evolution variables) is lost after that time. This 
may be due to reflections from mesh-refinement bound- 
aries. However, for both the <j> and x runs the errors are 
extremely small, and of comparable magnitude. 



FIG. 6: Errors in the real part of the I = 2, m = mode of 
r^4, extracted at r = 30toi. The upper plot shows results 
from runs with the <j> method, while the lower plots shows 
results with the \ method. See text for the grid details and 
discussion. 



on runs that use the initial-data parameters of the run 
"Rl" in 0], f° r which comparison simulations were also 
performed in We evolve these data with both the 

<f> and x variants of the moving-puncture method, and in 
each case compare runs with rj = 1 and rj = 2, to deter- 
mine which aspects of the simulations are most strongly 
affected by different gauge choices. 

The main goal of this paper is to demonstrate the ac- 
curacy and efficiency of our code, and of course to verify 
that it gives correct results. We begin by determining the 
grid setups necessary to achieve fourth-order-accurate re- 
sults, and present our results with error bars calculated 
from the difference between the highest-resolution runs 
and Richardson extrapolated values. By "grid setup" we 
do not simply mean "resolution" ; the sizes of the mesh- 
refinement boxes are also important, both for the accu- 
racy of the simulation, and the extracted physical quanti- 
ties. Having done that, we compare our extracted wave- 
forms with those produced by the independent LEAN 
code 0]. This is an extremely strong test: it validates 
both codes, and also demonstrates the high accuracy of 
their results. Finally, we study in detail the accuracy 
of various extracted physical quantities (radiated energy, 
angular momentum, and angular frequency during inspi- 
ral), and their dependence on the radiation extraction 
radius and the gauge parameter rj. 



VI. NUMERICAL EXPERIMENTS FOR TWO 
ORBITING PUNCTURES 

In this section we calibrate our code for binary evolu- 
tions. These will be the principal application of our code 
in future work, and we therefore perform a more detailed 
study than in the case of single black holes. We focus 



A. Setup 

The initial-data parameters for the runs in this sec- 
tion are: the punctures have mass parameters mi = 
■ni2 = 0.483, and are placed on the y-axis at y = ±3.257 
with momenta p x = =p0.133. The individual black-hole 
masses, as determined by Eq. |J5}, are Mi = M2 = 0.505, 



13 



and the total ADM mass of the spacetime is Eadm — 
0.996. These parameters correspond to the run "Rl" 
in 10]. They result in ~ 1.8 orbits before merger at 
roughly 160Af . We define three times indicating the 
merger time: tAH, the time when an apparent horizon 
first forms, t a , the time at which the lapse at the center 
drops below the value 0.3 (following e.g., 0>Is3)j an< ^ 
tmax, the time at which |^4| reaches a maximum (which 
depends on the extraction radius r ext ). While tAH is of 
immediate relevance regarding the simulation, it is also 
more costly to evaluate accurately, while accurate eval- 
uation is trivial for t a and t max . We therefore find it 
very useful to check convergence in phase by evaluating 
t a and t max , and note that t a and (t max - r ext ) give an 
estimate for tAH which is accurate to a few M. 

Note that in this section, all distances and times are 
either scaled with respect to the total black-hole mass, 
M (consistent with our discussion of single black holes 
in Sec. [Vj, in which case the appropriate unit is given 
(e.g., r^i{M or, when no rescaling has been done, 
the numerical coordinates are used (e.g., "extraction at 
r = 30"). 

We label the grid setups for orbit runs with the no- 
tation X[ni x Ni : n 2 x N 2 : buj][h~^ in : h max ], where 
X denotes the choice of conformal factor 4> or x, and 
the grid is composed of ri\ levels of Nf grid points and 
H2 levels of N$ grid points (reducing the number of grid 
points appropriately when discrete symmetries are ap- 
plied), and buf mesh refinement buffer points are used 
(not counting ghost zones for parallelization) . The quan- 
tities h m in and h max denote the grid spacing on the 
finest and coarsest levels. The qualifier X occasion- 
ally carries subscripts specifying further parameters. Ex- 
amples would be 0[5 x 32 : 5 x 64 : 6] [38.4 : 8] or 
X*7=o.05[5 x 32 : 5 x 64 : 6] [38.4 : 8]. The ratio of grid 
spacings between neighboring levels is always two. 

We have performed a large number of runs, both com- 
plete convergence series and lower resolution "exploration 
runs" with different grid layouts, gauges or numerical 
methods — for the presentation here we have to make 
a selection and present results from three series of runs, 
which we have found typical: 

1. BAM01: <j> n= i[4 x i : 6 x 2i : 6]. 

2. BAMxl: x„=i[5 x i : 4 x 2i : 6]. 

3. BAM%2: Xri=2[5 X i : 4 X 2i : 6], i.e., as above, but 
with rj = 2. 

The BAM01 series is representative of our early experi- 
ments. Apart from using the <f> evolution variable, they 
also used the ttt gauge advection choice. We later found 
the BAMx runs to be more accurate. In addition, the 
merger times converged from below, and convergence be- 
havior was monotonic even at low resolutions. For the 
runs presented here, we see no strong difference between 
the ttt and 000 gauge advection choices as described in 
Section Hi Dl However, the manifestly strong hyperbolic- 
ity of the 000 gauge j3|J makes that choice more attrac- 



tive. For the ttt choice, the slow-speed modes described 
in |35l ] can be clearly seen in animations of the grid vari- 
ables. Both choices yield stable evolutions, but we regard 
the 000 choice as superior and have used it in the BAMx 
runs presented here and in subsequent work. 

The runs presented here have 9 and 10 refinement 
levels (labeled from for the coarsest to 8 or 9), and 
use twice the number of grid points on the outer levels 
than on the finest levels, as detailed in Tabled We find 
that this setup yields higher accuracy for wave extraction 
without too drastic an increase in computational cost. 
Typical performance numbers of our code are displayed 
in Table [n] All runs are carried out with the symmetry 
(x, y, z) — > (— x, — y, z) and (x,y,z) — ► (x, y, — z), reduc- 
ing the computational cost by a factor of four compared 
with runs that do not exploit any discrete symmetries. 
The courant factor C — At /hi is kept constant, and is 
set to C = 1/2 for the inner grids, while for the outer 
grids at levels 0-2 the time step is kept constant at the 
value of level 3. All runs presented here use the RK4 time 
integration scheme. Using the ICN scheme (without ar- 
tificial dissipation) did not change results significantly. 
We find that for a constant courant factor we occasion- 
ally encounter numerical instabilities in the outer regions 
of the simulation domain, but these were cured by freez- 
ing the size of the time step in the outermost 3 (BAMx 
runs) or 4 (BAM(/>1 run) levels. 

All the BAM runs presented here use six AMR buffer 
points (see Sec. I1V|) . which is less than required to iso- 
late the fine level "half" timestep from time interpola- 
tion errors at the mesh-refinement boundary, and in par- 
ticular also less than required for the fully fourth-order 
Christmas-tree scheme suggested in [6^ . We have exper- 
imented with using higher numbers of buffer points up to 
the number required for the Christmas-tree scheme, but 
have not found significant improvements in the results, 
which is consistent with the fact that we find fourth- 
order convergence and no significant improvement of the 
results when decreasing the timestep. We conclude that 
at the resolutions presented here, six buffer points are 
enough to suppress errors from interpolation in time at 
mesh-refinement boundaries below the relevant thresh- 
old as far as the dynamics and low frequency waves 
are concerned. To suppress high-frequency reflections at 
the mcsh-rcfincmcnt boundaries, which can we have seen 
in quantities like ^4 or the constraints, we use fourth- 
order Kreiss-Oliger dissipation as described in Section llVl 
where the factor a is chosen as a = 0.1 in the inner lev- 
els and a = 0.5 in the outer levels (where the waves are 
extracted) . 



B. Results 

We have obtained fourth-order convergence for r^4 
and the puncture tracks for sufficiently high resolution 
in the BAM</>1, BAM^l and BAM^2 series, requiring at 
least 48 grid points on the fine levels. For the BAM(/>1 
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Run 




' ''max 


fmax 


0„=i[4 x 32 : 6 x 64 : 6] 


1/25.6 


20 


648.0 


<j> v= i[4 x 40 : 6 x 80 : 6] 


1/32.0 


16 


672.0 


<j>r,=i [4 x 48 : 6 x 96 : 6] 


1/38.4 


40/3 


666.7 


<^ =1 [4 x 64 : 6 x 128 : 6] 


1/51.2 


10 


680.0 


(/>,,=! [4 x 72 : 6 x 144 : 6] 


1/57.6 


80/9 


644.4 


X„=i, 2 [5 x 32 : 4 x 64 : 6] 


1/25.6 


10 


325.0 


X„=i,2[5 x 40 : 4 x 80 : 6] 


1/32.0 


8 


324.0 


X„=i,a[5 x 48 : 4 x 96 : 6] 


1/38.4 


20/3 


323.3 


X„=i, 2 [5 x 56 : 4 x 112 : 6] 


1/44.8 


40/7 


322.8 


X„=i,a[5 x 64 : 4 x 128 : 6] 


1/51.2 


5 


322.5 


X„=i,2[5 x 72 : 4 x 144 : 6] 


1/57.6 


40/9 


322.2 



TABLE I: Grid setups used for binary evolutions. See text 
for definition of the notation in the "Run" column. h m in 
and hmax (rounded to 3 digits) denote the finest and coarsest 
grid spacings, and r max is the location of the outer boundary 
(rounded to 4 digits). 



grid configuration 


procs. 


mem. 


(GByte) 


time 


M/hour 


X[5 x 56 : 4 x 112 


6] 


10 




8.9 


192 


18.2 


X[5 x 64 : 4 x 128 


6] 


12 




11.8 


306 


13.7 


X[5 x 72 : 4 x 144 


6] 


14 




17.5 


505 


9.7 



TABLE II: Typical performance results for runs lasting 
350Af : number of processors, maximal memory requirement 
in GByte (to be precise, we quote the resident size of the pro- 
gram, i.e. the physical memory a task has used), total runtime 
in CPU hours and average speed in M/hour for the Altix 4700 
of LRZ Munich [jjj (using Intel Itanium2 Madison 9M CPUs 
running at 1.6 GHz). 

series, low resolutions, e.g., with (f> v= i [4 x 32 : 6 x 64 : 6], 
i = 32,40,48, show no systematic convergence behavior, 
while for the x series, low resolutions show a convergence 
behavior as illustrated in the lower panel of Fig. [7| con- 
vergence is at least monotonic, but only the runs with 
high resolution are consistent with fourth-order conver- 
gence. 

We now focus on the BAM%2 series. The runs x»?=2 [5 x 
i : 4 x 2i : 6] show approximately second-order con- 
vergence for the set i — 32,40,48, but achieve clean 
fourth-order convergence with % — 48, 56, 72. The con- 
vergence results are shown in Figures [7\ and [S] Note 
that we see convergence in the waveforms without the 
time shift performed in |lfjj |. In particular, the bottom 
plot in Figure [7| shows the convergence of the time of 
the maximum amplitude in r^4 (extracted at 40M), as 
a function of the resolution on the finest box. We see 
that t max = 204.65 ± 1.3M for the second-order-accurate 
results, and t rnax = 203.9 ± 0.2M for the fourth-order 
accurate results. These results are extre mely accurate; 
compare, for example, with the results in |lOllll| . 



A few comments about these results are in order. 
One might have stopped at the second-order convergent 
i = 32, 40, 48 results indicated in the lower panel of Fig. 
[7| for the BAMx2-runs, since the code does have second- 
order ingredients (the initial data calculation and inter- 
polation in time at mesh-refinement boundaries). How- 
ever, this theory is vetoed by finding that neither the 
accuracy of the initial data, nor the number of mesh- 
refinement buffer points, nor the size of the time step, 
have a significant influence on the result. In such a sit- 
uation it is necessary to increase the resolution in con- 
vergence tests, and indeed for higher resolutions fourth- 
order convergence was found. Note also that fourth-order 
convergence does not extend to the early (before approx- 
imately t = 125M) and very late (after approximately 
t = 310AT) parts of the waveform, where the errors are 
very small (Fig. UJ, and to the late part of the puncture 
tracks after the merger as shown in Figs. 1161 1171 and ITH1 
The late-time loss of convergence may be due to the lo- 
cation of the outer boundary, which is at ~ 320M, and 
we expect incoming errors to reach the extraction radius 
(at 40M) after about 280M. It is encouraging that these 
errors are so small that they are detectable only in a con- 
vergence test, and that they do not have any effect on the 
stability of the runs. 

Having established fourth-order convergence of our 
runs in the given regime, and having determined a grid 
and parameter setup that produces accurate results, we 
now perform an independent validation of our results by 
comparing the ^4 waveforms with those obtained with 
the LEAN code, as published in I n the notation 

introduced above, the grid specifications of the high- 
resolution LEAN code run are ^=i[2 x 72 : 6 x 130 : 
3] [32.92 : 3.47]. This comparison is shown in Figs. and 
ITU1 Figure EH shows highest-resolution BAM and LEAN 
r\&4,? = 2, to = 2 modes, extracted at r = 30, for full 
runs. Figure focuses on the main part of the wave- 
form, and error bars are shown for the BAM results. We 
see that the results from the two codes show excellent 
agreement. 

The main focus of this paper is to present and validate 
our code, and show that we are able to produce highly 
accurate results with moderate computational resources. 
In this spirit we present numerical error bars for our re- 
sults, which are easily determined from the difference 
between highest resolution runs and Richardson extrap- 
olated values. However, waveforms and radiated energies 
also come with errors due to the finite extraction radius 
of the waves, and due to physically inappropriate bound- 
ary data. Fig. llll shows radiated energies versus time for 
X n =2 BAM runs Xj?=2[5 x i : 4 x 2i : 6], (i = 56, 72) and 
extraction radii r = 25,30,35,40. The dashed lines are 
from the i = 56 run and the full lines from the i = 72 run. 
The results have been shifted in time by the differences of 
extraction radii to minimize the phase difference. Clearly, 
the error from the variation of extraction radius is larger 
than numerical error at radii less than 30M. Assuming 
that the error falls off with some power of r, a curve fit 
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FIG. 7: Convergence of the BAM Xv=2[5 x i : 4 x 2i : 6][i = 
32, 40, 48, 56, 72] series. Top image: scaled differences between 
results for different resolutions demonstrate fourth-order con- 
vergence. Bottom image: the time t max (in units of M) of 
the maximal amplitude in 'M is plotted versus resolution, and 
Richardson extrapolated values are shown. The three most 
accurate runs (i = 48, 56 and 72) show fourth-order con- 
vergence, while the runs with (i — 32, 40, 48 and 56) show 
second-order convergence. 



of our results suggests that the error falls off as r 2 , and 
that in the r — > oo limit the extracted energy is 3.52%. 
The value extracted at r — 40 therefore has an error of 
only about 2%, in contrast to the numerical error, which 
is less than 0.3%. Further progress in accuracy obviously 
makes it very desirable to better model the fall-off prop- 
erties of the radiation. Although the only completely 
aesthetically plea sing solution would be to compactify at 
null infinity (see [zEEl El El El El E3 for some re- 
cent work and overviews) , one should expect that simple 
estimates based on perturbations of Kerr can be used to 
obtain significant improvements. A more detailed analy- 
sis goes beyond the scope of this paper, which focuses on 
the numerics, and will be published elsewhere. 

The idea of extracting wave signals at finite distance 
from the source is that the timelike cylinder traced out 
by a sphere at sufficiently large distance from the source 
can be viewed as an approximation of null infinity, and 
increasing the distance will increase the quality of the 
approximation. Thus an approximation to the signal ex- 
pected at a detector located at an astronomical distance 
from the source can be calculated. When the distance 
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FIG. 8: Convergence of the BAM Xv=^[5 x i : 4 x 2i : 6][i = 
48, 56, 72] series. Upper image: the early part of the waveform 
does not show fourth-order convergence until approximately 
t — 125M. Lower image: clean fourth-order convergence is 
lost at approximately t = 310M. 
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FIG. 9: Overlay of i?e(r* 4 ), extracted at r = 30 for high- 
est resolution LEAN <fi n= i (as published in |l3j0 and BAM 
X))=2 [5 x 72 : 4 x 144 : 6] [44.8, 5.714] runs. The results have 
been shifted in time so max(|\I r 4|) is aligned with t = 0. 



is of cosmological scale, cosmological redshift or further 
effects would have to be added "by hand". This idea 
raises several serious issues: the error introduced by the 
cutoff at finite radius needs to be estimated, "extrapo- 
lation procedures" to larger radius may yield significant 
improvements if a fall-off law for the finite distance re- 



16 



0.08 
0.06 
0.04 

sr 0.02 

3 o 

-0.02 
-0.04 
-0.06 



BAM rich+error / 

BAM rich-error /' 

LEAN high res. / 

// 

// 


\\ 

V 

V / 

\ / 


v / 

\ A 

\ / 
v. y ' 


\ / " 
xj 



-10 





Time (M) 



10 



FIG. 10: The Richardson-extrapolated "main" part of the 
waveform is shown for the BAM x»;=2 runs, error bars are ob- 
tained from the difference between Richardson-extrapolated 
result and the highest resolution run. The BAM result is 
overlaid with the highest resolution LEAN result. The re- 
sults have been shifted in time so max(|^4j) is aligned with 
t = 0. 
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FIG. 11: Radiated energies plotted versus time for the runs 
^=2(5 x i : 4 x 2i : 6] (i = 56, 72) and extraction radii r — 
25, 30, 35, 40. The dashed lines are from the i = 56 run, the 
full lines from the i — 72 run. Larger extraction radii yield 
smaller values for the radiated energy. The left image shows 
the result for the complete run, the right images zooms in on 
the late stage of the run. Results have been shifted in time by 
Eq. Jgnj to minimize the phase difference. Clearly, the error 
from the variation of extraction radius is not smaller than 
numerical error. 



suits can be assumed, and the gauge dependence of the 
results at finite distance needs to be addressed in order 
to optimize such procedures. It is therefore very valu- 
able to compute results characterizing the asymptotics of 
the gravitational field by different methods, which may 
show different effects from gauge, or different fall-off laws, 
and compare the results of such different prescriptions of 
asymptotic quantities. 

Along these lines, a good check on the consistency of 
the wave extraction algorithm is to compare radiated en- 
ergies with the energy balance that can be determined 
from evaluating the mass integral Eq. (|54|l at the begin- 
ning and end of the integration time. At finite extraction 
radius one can determine an estimate for the difference in 
Bondi mass; see Section fill Bl We have done this for all 
of our runs, and find excellent agreement for r\ = 1, and 
less accurate results (approximately 4% of radiated en- 
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FIG. 12: Shape of the apparent horizon in the x-y plane 
plotted each 10M for a simulation with initial separation of 
r — 3.5M. The dotted line represents the trajectory of the 
puncture. The behavior of the second horizon can be ob- 
tained by the symmetry of the problem. The common appar- 
ent horizon (peanut-shaped in the figure) appears when the 
black holes have merged. 



ergy, i.e., roughly 10% error) for rj = 2. The poor results 
for 7/ — 2 may be due a drift in the coordinates. Such 
a drift was seen in the coordinate radius of the horizon 
in the Schwarzschild case in Section IV Al For the case of 
orbiting punctures, we locate apparent horizons using a 
horizon finder based on the AHFinder code in the Cactus 
infrastructure . Figures El an d El show the motion 
of the apparent horizons for an orbital evolution (with 
punctures evolved from initial positions y = ±3.5), and 
also the coordinate radius of the single, and eventually 
common, apparent horizons as a function of time. We 
once again see an 77-dependent coordinate drift, which we 
expect also affects the quality of the Bondi mass. How- 
ever, note that such a strong gauge dependence is not 
seen for "J 4. 

Important quantities to be determined are the amount 
of radiated angular momentum, and the final spin of the 
black hole. For the time development of the angular mo- 
mentum as determined from the surface integral Eq. 156|) 
see Fig. 1141 We determine the initial angular momen- 
tum by means of the surface integral Eq. (|56|l . which 
can be evaluated analytically for Bowen-York data as 
Jo = pD = 0.866, corresponding to a Kerr parameter 
of a/Mo = 0.87. We numerically calculate the final an- 
gular momentum as J final = 0.634, corresponding to a 
Kerr parameter of a/Mfi na i = 0.688 and radiated an- 
gular momentum of 25%. The final angular momentum 
has been estimated by several methods with an error of 
roughly 1%. The surface integral Eq. ijSlfl) gives results 
that are very accurate and consistent between different 
choices of the ^-parameter. We also examine the complex 
quasi-normal ringdown frequency. The real and imagi- 
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FIG. 13: Coordinate radius tah of the apparent horizons 
as function of time for rj = and r\ = 2. The choice of r) 
has a strong influence in this gauge dependent quantity. The 
apparent-horizon mass Mah does not show any such gauge 
dependence. 
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FIG. 14: Angular momentum computed at r — 30 for gauge 
parameters rj = 1, 2 with the BAMx series — no significant 
dependency on the gauge is seen. 

nary part of the quasi-normal ringdown frequency can be 
determined by directly fitting an exponentially damped 
sinusoid to the gravitational wave signal. The imaginary 
part (quality factor) can easily be read off from |\I>4|, 
which shows exponential fall-off, and the real part of the 
frequency can be determined from the time derivative of 
the phase angle as defined in Eq. (|53|l . All these results 
are consistent, and yield a final spin of the black hole of 
J = 0.683Mj ina[ ± 1%, where Mf ina i is the mass of the 
final black hole. Remarkably, the orbital frequency of the 
punctures levels off (see Fig. IT3|) to the real part of the 
quasinormal frequency at late times. Note that at early 
times, before the merger, the wave frequency has been 
observed to be twice the orbital frequency [H| , as would 
naively be expected from the quadrupole formula. 

One of the remarkable facts about recent simulations 
of binary black holes, whether done in a generalized har- 
monic or BSSN moving punctures framework, is the qual- 
ity of the coordinate conditions: not only do they produce 



a "nice" spiraling motion with almost spherical (apart 
from a short time during merger) apparent horizons as 
seen in Fig. ^1 but the coordinate tracks give rise to a 
good estimate of the waves via the quadrupole formula 
(see, for example, |79)'). and the measured angular veloc- 
ities coincide very accurately with what is expected on 
physical grounds. For example, at the beginning of the 
simulation the angular velocity quickly reaches a value 
close to that expected from the initial data (approxi- 
mately Mil = 0.05); see also 0. A heuristic explana- 
tion for the latter fact has been given in |33| : symmetry- 
seeking gauge conditions (e.g., the gamma- freezing con- 
dition used here) should be expected to find an approx- 
imate helical Killing vector. This Killing vector will be 
unique up to a rigid rotation of the form uj x r . We 
can choose either corotation, or vanishing rotation via 
the shift boundary condition at infinity. In this work 
the shift is set to zero at infinity. Thus the punctures' 
coordinate speeds are expected to be equal to their phys- 
ical speeds seen from infinity. It is interesting to check 
whether the choice of gauge - in our approach this boils 
down to the choice of shift damping parameter 77 — has 
an effect on the coordinate tracks of the punctures. We 
plot the radial and angular motion of the punctures with 
numerical errors (again determined from the difference 
of the Richardson extrapolated value to highest resolu- 
tion result) for the BAMx2 and BAM01 runs in Fig. US 
While the results look consistent between different runs, 
a small but significant ^-dependent deviation can be seen 
in the radial motion, whereas a difference in angular mo- 
tion is not visible on the plot. 
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FIG. 15: Orbital motion with numerical errors obtained from 
difference of Richardson extrapolated value to highest resolu- 
tion result (top left: r(t), top right: ui(t)). The dashed curves 
represent Xv=2 and the full curves cj>ri=i- Bottom: The log- 
arithm of the puncture radial position shows an exponential 
decay at late times, with a damping time of r = 3.48 ± 0.01 

(Mfinai). 

Fourth-order convergence of the puncture motion in 
the %2 and cj)l runs is demonstrated in Figs. ^0 d an d 
ITH1 The coordinate angular speed u> is calculated using 
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u> = \/3\i/r, where \{3\i is the norm of the shift vector eval- 
uated at the ith puncture that gives the coordinate speed 
of the puncture across the grid. This quantity shows 
fourth-order convergence up to the black-hole merger. 
After that time the puncture distance from the origin 
decays exponentially (see Fig. I15J) . thus rather soon the 
punctures are both less than one grid point from the ori- 
gin and the convergence in oj deteriorates to first-order. 
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FIG. 16: Fourth order convergence for r(t) demonstrated for 
X>7=2 and (j>ri=i series (left and right upper images). Lower 
images: relative errors obtained from difference of Richardson 
extrapolated value to highest resolution result (left: x>7=2, 
right: (t> v =i). 
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FIG. 17: Fourth order convergence for demonstrated for 
X>7=2 series, at late times only first order convergence is seen. 
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FIG. 18: Fourth order convergence for ui{i) demonstrated for 
4>n=i series, at late times only first order convergence is seen. 



VII. QUASI-CIRCULAR ORBIT PARAMETERS 
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FIG. 19: Relative errors obtained from difference of Richard- 
son extrapolated value to highest resolution result (left: x»7=2, 
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tion runs. These are in turn based on those from Cook's 
1994 initial-data study [8(j- In that work the parameters 
for quasi-circular orbits were determined using an "ef- 
fective potential" (EP) method, whereby quasi-circular 
orbits corresponding to a given total orbital angular mo- 
mentum J were identified by the minimum in a curve of 
the binding energy = Eadm — Mi — M 2 versus the 
proper separation, in analogy with Newtonian physics. 
Cook's results applied to inversion-symmetric (not punc- 
ture) Bowen-York data, but a later study of the inner- 
most stable circular orbit (ISCO) of puncture Bowen- 
York data by Baumgarte [81| suggested that there are 
only minor (if any) physical differences between the two 
types of initial-data set. 

Note that Cook's parameters cannot be directly used 
for punctures. For this reason, the parameters we actu- 
ally use come from Baker et al. (8?J who have translated 
Cook's parameters into puncture parameters. As pointed 
out in |84j , this translation has introduced additional er- 
rors on the order of 1% in the masses. The translated 
Cook parameters have been found in practice to produce 
reasonably convincing quasi-circular orbits. (However, 
see PJJ, |23 for evidence and estimates of eccentricity for 
larger initial black-hole separations.) In addition, Baker 
et al. found that the merger waveform was largely inde- 
pendent of the initial separation of the punctures, sug- 
gesting that the parameters predicted by the effective- 
potential method really do correspond to points on an 
inspiral sequence. 

However, there are a number of alternative ways to es- 
timate the momenta as a function of separation for black 
holes on an inspiral sequence. One option is to use an 
approach based on the assumption of the existence of 
a helical Killing vector (HKV) miHH, for which a 
sequence of parameters for puncture data has been com- 
puted |2?1 | . Another option is to use parameters predicted 
from post-Newtonian (PN) theory. How sensitive is the 
final waveform to each of these approaches? 

Let us first consider post-Newtonian parameters. For 
a given separation D, the momentum of each puncture 
can be given to 3PN order in the ADMTT gauge by H3 



For the runs in Section IVII we chose the same initial 
parameters as used by Baker, et al. |lfj |. for their calibra- 
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Method 


m/M 


y/M 


p/M 


Effective-potential (EP) 


0.4782 


3.2248 


0.1317 


Helical Killing vector (HKV) 


0.4782 


3.2248 


0.1307 


3PN 


0.4780 


3.2248 


0.1329 



TABLE III: Initial parameters for a given initial coordinate 
separation, from three different approaches: the effective- 
potential (EP) and helical Killing vector (HKV) methods, and 
the 3PN formula 1641 1. The parameters are scaled with respect 
to the total black-hole mass in the initial data. 



(163tt 2 - 4556) v + 104i/ 2 ] ( — 



~D 



7/2 



(64) 



The total mass is M — Mi + M2, the reduced mass is 
\i = MiM 2 /M, v = (J./M, and the PN order of each 
term is indicated by e. For equal-mass black holes with 
Mi = M 2 = 0.5, we have fi = v = 0.25. 

Eq. (|fj4"|) was derived using equations (5.1)-(5.3) in [85|. 
and noting that uj sta tic = and ^kinetic = 41/24 |8g. 
Equations (5.1) and (5.3) can be rearranged, order by or- 
der in the post-Newtonian expansion, to give the orbital 
angular momentum J as a function of puncture separa- 
tion D, and we then use the relation J = pD (which 
holds by definition; see section IV of jg^) to write p as 
a function of D. Eq. Ij64(l is not gauge invariant, but the 
ADMTT gauge is expected to be very close to the con- 
formally flat gauge that we choose for our initial data. 

Initial parameters from these three approaches, 
effective-potential (EP), helical Killing vector (HKV), 
and third-order post-Newtonian (3PN), are given in Ta- 
ble IIIII The parameters are scaled with respect to the 
total black-hole mass, M = M\ + M2, which is a con- 
venient quantity in all three approaches; in the post- 
Newtonian expression Q64JI the black-hole masses appear, 
but the other standard mass scale, the total ADM mass 
Madm-i does not. In each case, the coordinate separation 
of the punctures is kept fixed, and a prediction for the 
momenta that will produce a quasi-circular orbit is pro- 
vided. These predictions all differ by less than 1% (which 
is also the error estimate in Cook's sequence an d 
we might expect the resulting orbital motion and merger 
waveforms to be equally close. However, from Figure |2*UI 
it is clear that this is not the case. 

The black-hole merger times differ by about 40M, and 
the merger waveforms arc noticeably different. However, 
from the tracks of the puncture locations during the evo- 
lutions, it is not clear which is the "better" choice of 
quasi-circular orbit parameters, or even what it would 
mean for one choice to be better than the other. These 
evolutions also suggest that parameters calculated from 
PN methods are an acceptable alternative to numerically 
generated parameters, allowing a wide range of configura- 
tions to be explored without the need for accompanying 
initial-data studies. 
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FIG. 20: Results from evolution of Bowen-York puncture data 
using three choices of initial parameters, described in the text. 
Top: real part of the I = 2,m = 2 mode of r^i, extracted at 
r = 40M. Lower left: paths followed by the punctures during 
evolution. Lower right: energy extracted at r — 40M after 
300 A/ of evolution. 
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FIG. 21: Results from evolution of Bowen-York puncture data 
using three choices of initial parameters, described in the text. 
Left: Absolute value of the I = 2,m = 2 mode of r^4, ex- 
tracted at r = 40 M. Right: Phase angle (in units of 2n, 
shifted in time to align the maxima of afos^^) at t = 0, and 
also aligning the phase at t — 0. 



Despite the apparent differences in the dynamics, the 
radiated energy from the merger differs by only a few 
percent. In Figure I2T1 the amplitude and phase of ^4 are 
shown separately, with the results shifted in time so that 
the maximum in the amplitude occurs at the same time 
for all three choices of initial parameters (see [lbll87|| L It 
is now clear that, although the dynamics and waveforms 
look quite different, these differences are merely cosmetic: 
the physics of the merger is the same for all three choices 
of initial parameters. 



VIII. DISCUSSION 

We have presented a new code to evolve black-hole 
binaries, which is an extension of the older BAM code 
[6L ll5L fl6j | . The new BAM code implements the "moving 
puncture" method (in both its <fi and X versions) within a 
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moving-box-based adaptive mesh-refinement grid struc- 
ture, and uses fourth-order-accurate spatial finite differ- 
encing and RK4 time evolution. The primary analysis 
tool, the extraction of gravitational radiation waveforms, 
is implemented using the '5 4 Newman-Penrose scalar. 

In this paper we have presented a number of impor- 
tant tests of our code: evolutions of a single non-spinning 
black hole, a single spinning black hole, and, the cru- 
cial test, evolutions of black-hole binaries. In addition 
to demonstrating fourth-order convergence in regimes of 
physical interest, each test has provided valuable insight 
into the grid sizes, resolutions and geometries necessary 
to achieve accurate and efficient simulations. One of the 
ultimate uses of our code will be to perform large param- 
eter studies of gravitational- wave sources, and therefore 
efficiency of the code is of equal concern to its accuracy 
and stability. The performance of the code, and how to 
tune its configuration to obtain good performance, is one 
of our main results: as shown in Table ITT1 we are able to 
perform a convergence series of three runs in the fourth- 
order convergent regime at a total cost of roughly 1000 
CPU hours and a maximal physical memory consump- 
tion of less than 18 GByte. Accurate results for binary 
black hole evolutions can thus already be obtained with 
relatively small commodity clusters, which are available 
in many research groups. 

Evolutions of single black holes showed that both the 4> 
and x moving-puncture methods are stable and accurate, 
although the \ method shows better convergence proper- 
ties at the puncture. The moving-puncture method was 
also found to be stable even when resolutions of up to 
M/512 were used at the puncture. For black-hole bina- 
ries, the x method showed monotonic convergence behav- 
ior (for example, in the merger times) at low resolutions, 
and appears to us as clearly preferable, even though we 
have not found prohibitive problems with the <j) method. 
As a further element of validating our code, we have com- 
pared our waveform data with that from the independent 
LEAN code , and found that the two were in excellent 
agreement. This is a strong validation of both codes, and 
it will be interesting to see more extended comparisons 
between codes capable of performing binary black hole 
simulations - both regarding their results and efficiency 
in order to further refine the methods of the field. 



We have investigated the influence of gauge choice, and 
have found consistent results between the 000 and ttt 
shift advection choices, as well as the choice of the r\ pa- 
rameter in the T-driver shift condition, e.g., also by com- 
paring the LEAN and BAM results. We have also found 
that, as expected, the rj parameter in the F-driver shift 
condition has an effect on the final coordinates of the so- 
lution, and as a result larger values of rj lead to a larger 
coordinate size of the black hole; in effect, the black hole 
becomes better resolved on the numerical grid. This co- 
ordinate drift may however also cause problems for naive 
wave extraction algorithms, and further research will be 
required to make optimal gauge choices. 

Finally, we performed simulations using different 
choices of initial parameters for the momenta of the punc- 
tures. We found that very small changes in the initial 
momenta can make a large difference in the merger time 
of the black holes, but do not change the physical prop- 
erties of the radiation. 

Having carefully tested and calibrated our code for 
simulations of comparable-mass black- hole binaries, in 
future publications we plan to extend our research to pa- 
rameter studies of unequal mass and spinning black holes, 
and to use initial-data sets with larger separations. 
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